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Abstract. We introduce a nonlinear modification of the classical Hawkes process, which allows inhibitory cou- 
plings between units without restrictions. The resulting system of interacting point processes provides a useful 
mathematical model for recurrent networks of spiking neurons with exponential transfer functions. The expected 
rates of all neurons in the network are approximated by a first-order differential system. We study the stability 
of the solutions of this equation, and use the new formalism to implement a robust winner-takes-all network that 
operates robustly for a wide range of parameters. Finally, we discuss relations with the generalised linear model 
that is widely used for the analysis of spike trains. 

1. Introduction 

The problem of formulating and investigating mutually interacting point processes is of great importance 
both in the theory of point processes and in their applications. The classical model is due to Hawkes [9l [10] . He 
considers a point process that is defined by a specification of its rate function X(t) (called "intensity" in Hawkes' 
papers, or "conditional intensity" elsewhere in the mathematical literature). The value X(t)St is the expected 
number of events in the interval (t, t + St), and the rate itself is defined as a random variable obeying the dynamic 
law 

A(t) :=A 

tj<t 

where tj is the time-stamp of the jth event, and K is a, positive kernel to ensure positive rate. Of course 
it is possible to choose X(t) as a vector, and to allow its components to be influenced by events in the other 
components, assuming a separate kernel for each pair of components. In this way, one obtains a family of 
processes that interact linearly. Applications of Hawkes' theory in seismology have been quite successful, see [17] 
for a review. Applications in the neurosciences, however, are rare, see |12j and references therein. This is mainly 
due to the fact that positive kernels only allow one to model mutual excitation. A fundamental feature of most 
biological neural networks, however, is the presence of inhibitory couplings. So, Hawkes' model falls short as a 
model for biological neural networks as it cannot represent retarding interactions. 

Here we propose an alternative model which goes beyond Hawkes' linear formalism, adhering to a representa- 
tion in terms of rates and avoiding to invoke secondary state variables like the membrane potential. Specifically, 
the change in the instantaneous rate due to an incoming event at time t is given by 

X(t + e) = wX(t), 

where w is the "weight" of the connection under consideration. In this framework, w > 1 yields an excitatory 
connection, w < 1 gives an inhibitory connection, and w aa ' — 1 means that the corresponding link is inactive or 
absent. Based on this principle, one can construct networks of computational units, each of them characterized 
by its own instantaneous rate A a (t) . The weights of all connections are encoded in a matrix (w aa ' ) of positive 
numbers. In the first part of the paper we are mainly concerned with the expected instantaneous rates EA a (t) = 
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y a (t), sometimes also plainly called "intensity" or "rate" in the literature. For the expected rates we are able to 
heuristically infer an ordinary differential equation that approximates the expected instantaneous rates 

dy a (t) Mi 
df = Va{t) y a > (t) log w aa > , 

a' 

and we explore its range of validity with numerical simulations. Similar models have been introduced in [23, 24 , 
based on a slightly different approach. Both approaches, however, have many common features with the class 



of generalized linear models introduced in [25] and [T5], and with a class of cascade models [18]; see Section 3.4 
for details. The similarities trace back to the fact that the multiplicative rule is additive in the logarithm of the 
instantaneous rates, which are a natural parameter (likelihood) for certain point process models. 

The description we choose is based on the inhomogeneous Poisson process, viewed as a continuous-time 
Bernoulli process. This is possible since the rate function X(t), i.e. the (normalized) expected number of events 
in the time interval (t,t-\-St), and the probability p(t) that the interval (t,t + St) contains at least one event, are 
connected by the relation 

p(t) = 1 - exp(-\(t)5t). 

If St is infinitesimally small, we have p(t) = X(t)St and it is possible to use the above expression to compute the 
expected value of the rate function. We decided to model the point process as a binary process on an infinitesimal 
grid. This approach is equivalent to the measure theoretic one by means of non-standard analysis, see the 
axiomatic treatment |15j . This approach has some advantages though: First, it is intuitive, mathematically 
rigorous and avoids measure-theoretic complications. Second, the non-standard infinitesimal discretization step 
used to derive theoretical results can alternatively be fixed as a small standard number, which in a natural way 
leads to a Monte Carlo simulation scheme. 

Our paper is organized as follows: In Section [2] we introduce the Cox process (doubly stochastic Poisson 
process) on an infinitesimal grid and establish some preliminary results. We further define multiplicatively inter- 
acting point processes and derive an approximate differential expression for the expected rates. In equilibrium, it 
corresponds to a system of ODEs that we call the rate equation of the system. In Section[3]we study transmission 
properties of single neurons, i.e. of "networks" consisting of one single Poisson input and an integrator. Further, 
we explain how our model is related to other common models in computational neuroscience. In Section [4] we 
investigate the rate equation of the system more thoroughly and systematically analyse small networks consisting 
of 2 units driven by Poisson input. We also show how it is possible to implement an efficient winner-takes-all 
dynamics in this framework. Finally, we discuss the scope of our results and indicate possible directions for 
future research in Section [5] 

2. Definition of the process and first-order properties 

Monte Carlo type simulations are of great importance in the study of stochastic processes, and they are usually 
performed on a discrete grid 

Mst :={kSt:ke N}, 

of resolution St, where St is a small positive number. On a mathematical level, this approach has the advantage 
that many results can be obtained by algebraic calculations. Then, the parameter St is sent to and, after 
verifying convergence conditions, the results can be transferred to the continuous-time stochastic process. 

One method to overcome certain technical issues and measure-theoretic complications when going to the limit 
of continuous time is to work on a grid 

H e := {ke:ke N*}, 

where now e is some infinitesimal number, and N* is the set of non-standard natural numbers, as in Nelson's 
internal set theory [15 . We will follow this approach to define interacting point processes, suppressing the explicit 
reference to e whenever possible. 

Now and in the rest of the paper, the reader not interested in the details of non-standard analysis should 
simply think of e as a really small number. As a matter of fact, all simulations were realized with such a scheme. 
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We refer to [2j EJ |T5j [16] for short introductions to the subject, and for a description of methods of non-standard 
analysis in the theory of stochastic processes. All tools of calculus we need in the paper are contained in [T5] , 

2.1. Cox processes on the grid. As a warm-up, and to fix some preliminary results we will need in the 
following, we define the Cox process and list some elementary properties of a Bernoulli variable driven by an 
infinitesimal positive random variable. For any two positive numbers x, y, we will use the notation x ~ y for 
expressing the fact that ~ is infinitesimal. 

Proposition 2.1. Let r be a positive random variable and X an independent Bernoulli random variable with 
parameter p = 1 — exp(— re). Then 

(2.1) EX ~ eEr, 
and also 

(2.2) E(l - exp(-re)) ~ eEr. 
Finally 

(2.3) Var(X) = EX(1 - EX). 

The proof of these facts is purely algebraic and can be found in the Appendix. We now move to the definition 
of a Cox process on the infinitesimal grid. To begin with, we recall that a grid stochastic process is a set of 
random variables (A(t))tgm indexed over the infinitesimal grid H. 

Given a positive grid stochastic process (A(t))tgn, a grid Cox process (X(t)) t ^m is an independent family of 
Bernoulli random variables, indexed over H, with time-dependent parameter 

p t = 1 - exp(-A(t)e). 

Finally, if there is a deterministic function fj,(t) on the infinitesimal grid such that X(t) ~ almost surely, then 
we call (X(t))t e n an inhomogeneous Poisson process. 

It is easily seen that this definition is equivalent to the standard definition of a Cox process. For instance, 
the random variables X(t) are independent Bernoulli variables, conditionally on their rate. We will prove that 
the expected count equals the integral of the expected rate. During the rest of the paper, the symbol (X(i))tea 
will denote a Cox process with rate A(i). In fact, the symbol A(i) denotes a positive stochastic process. For 
the Poisson process, it is possible to express the expected number of events as the integral of the rate function. 



Equation (2.1l yields 

EN(t) = E ^ X(s) = J2 Ex ( s ) - Y eEA ( s ) - / E H s )ds. 
se[o,t]n «e[o,t] H »e[o,t] H 

This proves 

Proposition 2.2. Denote by (N\(t)) t £u the counting process defined by 



N(t):= >^ X(s). 

Then 



E 

«e[o,t] H 



EN(t) ~ / X(s)ds. 
Jo 



The function N is not differentiable, so it does not make sense to consider the derivative =jrr. We introduce 
an operator ^ that acts on functions defined on EL 

Definition 2.3. If f : H — > R* is a function defined on the infinitesimal grid H £; then 

Af _ f(t + e)-f(t) 
At- e ' teM - 
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Of course, if f(t) is differentiablc in the standard sense, then 

dt At 

as has been proven in |15j . The following result will be used in later sections 
Proposition 2.4. The grid differential of the count process satisfies 

( 2 4) AN{t) ~ X{t + e) 

1 ' ' At e 

2.2. Multiplicatively interacting processes. We are now going to introduce a family of Cox processes which 
interact with each other on the basis of their events. To see how it works assume that X = (X a (t)) is a family 
of conditionally independent Bernoulli random variables with rates X a (t), indexed by some set A. In fact, even 
if the rates X a (t) are defined in terms of the realizations of X at times before t, the property 

EX a (t)X a ,{t') =EX a (t)EX a ,tf) = (l-exp(-A Q (t)e))(l-exp(-A ,(t')e)) = e 2 X a {t)X a ,{t') 

still holds, conditionally on the rates. 

Definition 2.5. Consider a positive coupling matrix W := {w aa >) and define rate functions by the relation 

(2.5) X a (t) := A a (0)exp( 

N a ,(t- e)\og 

Waa'j- 

a'eA 

The family X of the corresponding Cox processes are called multiplicatively interacting point processes with 
coupling matrix W . 

The stochastic time evolution of such process is captured by the random variables X a (t + e) which, for any 
time t and any given X a (t), satisfy the relation 

(2.6) X a (t + e) = A a (t)exp(^ X a ,(t) loguw)- 

a'eA 

We will refer to the variables A a as "instantaneous rates". In Hawkes' papers the same variables are called 
"intensities", whereas the expression "conditional intensities" is used in the mathematical literature. During the 
rest of the paper, the symbol X will denote a multiplicatively interacting family. 

2.3. Expectations. The aim of this section is to derive an approximate differential expression for the time- 
dependent expected rates. Once this expression is found, we experimentally show to which degree it predicts the 
equilibrium behavior of the stochastic system. The strategy is the following: 

(1) Derive an expression for the expectation of the grid differential, conditional on the actual rates. 

(2) Use this information to derive a differential expression. 

We stress that, conditional on the actual rates, the grid differential is a random variable which is independent 
of the actual realization of the point process, and which satisfies 

A> ^ | X a (t) ~ A 8 (t) EA a /(t)loguw. 

- 1 a'eA 

We are now almost in the position to derive the desired differential expression for the rates. Of course it cannot be 
expected that the expression we are looking for contains only expectations of rates; it turns out that covariances 
of pairs of rate variables also appear. For all a £ A the random variables X a (t) satisfy 

(2.8) AE ^ = J2 logw aa ,E{X a (t)X a ,(t)} ~ Y, ±ogw aa ,[EX a (t)EX a ,(t) + Cov(X a (t),X a ,(t))]. 

a'eA a'eA 



(2.7) E 



Assuming that the rates A a are (approximately) uncorrelated, Equation (2.8 1 can be used to guess a system of 
ODEs that describes the evolution of the event rates. In fact, we performed numerical simulations of networks of 
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FIGURE 1. Scheme of the stochastic Lapicque's integrator. 



sizes up to 100 neurons, both with specific architectures and with random topologies, specifically testing the cases 
which are critical for other type of processes. These simulations showed that the component processes indeed 
become uncorrelated after some time of relaxation, a finding which is supported by preliminary mathematical 
analysis involving covariances. As a consequence, we are convinced that the following definition is (heuristically) 
justified. 

Definition 2.6. (1) Define £ aa > := logw aa i . The system of ordinary differential equations 

(2-9) d ^=y a {t) YsVa'Waa', V(0) = Vo- 

a'eA 



is the rate equation associated with the system (2.5) 



(2) A family of interacting point processes is said to be in equilibrium if 

EX a (t) = y a (t), 

withy a (0)=E\ a (0). 

We stress that we of course have not proven that a family of interacting point processes always converges 
to equilibrium in the above sense. In fact, it is not even clear that interacting families in equilibrium exist at 
all. Again, extensive Monte Carlo experiments showed that EA a (i) indeed converges to the fixed point of the 
associated rate equation for large times t, and that interacting families indeed run into an equilibrium state after 
an initial transient. For the time being, a rigorous proof of this interesting numerical observation must remain 
open though. 

3. The stochastic perfect integrator 

Our goal is to study the behaviour of networks of multiplicatively interacting processes. Before we address 
this problem, we study the simple case of a single neuron which is fed with excitatory Poisson input. We call 
this very elementary system a stochastic perfect integrator, SPI in the following. 

Of course, the power of our model cannot be observed here, i.e. the possibility of modeling inhibitory synapses 
while keeping the mathematical analysis simple. However, it is useful to discuss this example to show what is 
the qualitative behaviour of the model, and to explore the connections with more established models. 



3.1. Adiabatic regime of the SPI. As we have already pointed out, Equation (2.9 1 does not predict exactly 
the behaviour of the rate dynamics. However, one could hope that, at the equilibrium, correlations do not play 
any role for the network dynamics. We call this regime as the adiabatic regime and we illustrate its features for 
a elementary system. Let us shortly illustrate its architecture. 
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Figure 2. Estimation of the instantaneous rate for a stochastic perfect integrator with different 
initial rates. Rate is estimated by convolution of spike data with a triangular kernel of width 
0.01. For the non-adiabatic simulation 10 5 trials were used and 5- 10 4 were used for the adiabatic 
simulation. Parameters are W22 — 0.01 and W21 = 1.2. Upper box: the initial rate is deter- 
ministically set to 1: the expectation E[A(t)|A2(0) = 1], < t < 2, Ai = 50 is plotted. Large 
oscillations due to the autocorrelation can be observed. Lower box: SPI with warm-up time. 
The expectation E[A(i)|A 2 (0) = 1], 15 < t < 17 is plotted, where Ai(i) = 25.26 for t < 15, and 
Ai(t) = 50 for 15 < t < 17. In this case the predicted firing rate yields a good approximation of 
the observed one. 



The system is composed of two units. The first unit has no self-inhibition, i.e. wn = 1, and it feeds input in 
to the second unit with a constant rate A and a weight W^x- The second unit has self- inhibition W22 and but no 
outgoing connection. 

Finally, we have to specify in which state we start the system. Let us first choose A2(0) = 1. The rate 
dynamics of the rate r is given by 

dr(t) 

- r{t){\ogW2i\ + \ogw 2 2r{t)) 1 



dt 

log t0 2 2 



and the right hand side equals if r(t) = — A |°^ ™^ . We define = ln(uiy ). It is a Riccati equation, the solution 
of which is given by 

A^ie A£2lt 



r(t) 



If A2(0) ^ 1, the equation can still be solved analytically. Now it is possible to compare the trajectories of the 
analytic solution with the trajectories of the expected firing rate in numerical simulations. It turns out that they 
do not coincide if the initial value of the rate is chosen to be exactly 1. Although, as shown in Figure [2j the 
observed average rate indeed converges to the fixed point of the rate equation, the precise orbit oscillates around 
the analytic solution. We stress that the initial value of the instantaneous firing rate of the output unit is fixed 
to A 2 (0) = 1, deterministically. The firing rate at equilibrium is — A }°|^ = 1.98. 
Summarizing, in Figure [2j upper panel, two different phenomena can be observed 

(1) the firing rate at the equilibrium is correctly predicted; 

(2) the transients oscillate around the analytic solution. 

We conclude that initializing the system on a given, deterministic value does not lead to a system in the adiabatic 
regime. 
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To solve this problem, let us observe that in the derivation of Equation (2.9 1, the variable 1/2 represents the 
expected value of the random variable encoding the rates. We conclude that we must choose the initial rate from 
the equilibrium distribution of the rates. 



Since the equilibrium distribution could not be obtained by analytic means, see also Section 3.2 we had 
to follow an alternative approach to obtain a reasonable solution. We describe in details the protocol of the 
simulation from which the plot in Figure [2] lower panel, was obtained: 

(1) we computed the input rate X wu such that the output rate at the equilibrium is 1 by the formula 

X W u = — , — " = 25.26; 

log W21 

(2) for 15 seconds we stimulated the output neuron with the rate X wu ; 

(3) at time we switched the input rate from X wu to 50. 



One now sees that the averaged spike histogram follows with very good accuracy the solutions of Equation (2.9 1, 
plotted in red. 

Finally, we want to spend some words about the following problem: Is it possible to map the parameter of 
the SPI to the parameters of a perfect integrator with Poisson input? A perfect integrator is characterized by 
a threshold T such that, if the membrane potential V raises above threshold, an output spike is emitted. We 
assume that each pre-synaptic spike produces an increase of the membrane potential of i. So, if the pre-synaptic 
spikes arrive with rate A, one sees that the output rate of the perfect integrator is given by A;p. Indeed, the 
stochastic perfect integrator has the same output rate as the deterministic perfect integrator if i = log W21 and 
T = — log W22 ■ The relation between the SPI and integrate-and-fire neurons is deeper than the pure possibility 
of mapping parameters of one model into parameters of the other: we will address this issue in more detail in 
Section 13.31 



3.2. Master equation of a stochastic perfect integrator. A full explanation for the observed transients 
can be given in terms of the evolution of the rate distribution. If the rate r has time-dependent distribution 
f(r,t), then the rate at time t is given by 



m = f 

Jo 



rf(r, t)dr. 



Deriving the master equation for the rate distribution is necessary to understand the system thoroughly (see 
Appendix [C| for a derivation). This equation reads 

(3.1) = —/(—,*) + — .*) - (* + r)/(r,t), 

at W21 W21 W22 W22 

complemented with the initial condition 

/(r,0) = /„(r). 

Here, A is the rate of the input process. A thorough analysis of this equation lies beyond the scope of this paper, 
but we would like to add some considerations. 

Let us first develop a heuristics for the asymptotic distribution of the rates. Assume that we initialize the 
system with a deterministic rate r(0), fix a small number 8t and denote by Ik, respectively Ok, the number of 
input, respectively output, spikes in the interval (kSt, (k + l)St]. After time t = KSt the rate will satisfy 

K 

(3.2) r(t) = r(0) ]1 ™ 2 V2 2 fe • 

fe=0 



Equation (3.2) shows that r(t) is the product of a sequence of random variables. Although these random 



variables are neither independent, nor identically distributed, one could hope that the logarithmic central limit 



theorem should hold in some weak sense. Of course, for deterministic initial rates Equation (3.2) shows that 



the distribution of r(t) will strongly oscillate, being only supported on a finite subset of R + , contradicting 
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FIGURE 3. Continuous blue line: estimated distribution of the rate variable at the equilibrium 
in semi-logarithmic plot. Dashed red line: Gaussian distribution fitted on the mean and the 
standard deviation of the final rates. During the simulation we kept track of the final values 
of the instantaneous rates after a long run (200s). The smoothed version was obtained by 
convolving 100 data points with a Gaussian kernel. 



the logarithmic central limit theorem. However, choosing r(0) from some distribution supported on the whole 
positive real line will avoid this effect. Simulations showed that the limiting distribution is a distorted lognormal 
distribution, indeed, in good accordance with the arguments we have just exposed. We have visualized the results 
in Figure [3] In the plot, a smoothed version of the empirical distribution of the logarithm of the final rates after a 
long run is shown. One sees that the distribution is a slightly distorted Gaussian, thus supporting our heuristics. 

Further qualitative evidence for the goodness of the lognormal approximation can be gained from the moment 
equation. This is derived in Appendix [D] and reads 

(3.3) = K2 - l]M«+i(i) - A[l - u&K(t). 



At the equilibrium we obtain the recursion 



w 22 1 



11 • 



which for large n behaves as 

H n+ i ~ Ait>2iMn = ^ ex P (nlog(w 2 i)) . 
On the other hand, the lognormal distribution satisfies 

/i n+ i = exp (2n + 1)J = exp ^"y ) exp ( na ' 2 ) ' 
From these equations we see that either distribution satisfies an asymptotic moments recursion given by 

Mn+i = fciexp(nfc 2 )^i„, 

with some positive constants k\,ki- This shows, that the tail scaling is very similar for the distribution of the 
asymptotic rates and for the lognormal distribution. 

We mention that these theoretical findings are in accordance with the fact that, in real neurons, the membrane 
potential is found to be normally distributed, see e.g. [5J. 
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3.3. Connection with leaky integrator models. We mentioned in Section [3~T| that the connection between 
the SPI and standard neural models goes beyond parameter mapping. In fact, Equation (2.9 1 can be obtained 
from the Lapicque's perfect integrator by the following method. Recall that a network of linear neurons can be 
described specifying the membrane potentials V a by the convolution 

V a (t)= X)(tfaa' **«')(*)• 
a'eA 

Here, the function K aa > descibe the post-synaptic potentials. Now, we assume that the neuron has transfer 
function F a , so that the instantaneous firing rate is given by X a (t) = F a (V a (t)) As a consequence, we obtain 



F'(V a (t))"^> = F'(V a (t)) £ (K' aa ,*X a ,)(t). 



dt y " v " dt 

a'eA 



For a perfect integrator, the kernel K aa i is the Heaviside function, and this has as derivative the Dirac S. The 
above equation then yields 

d\ a (t) 



F'{V a {t)) J2 «W (*(*)* *<»')(*)• 



dt 

a'eA 

We now choose an exponential transfer function F a (x) := exp(x) and obtain 

^l = Xa{t )J2^aa'X a ,(t). 

a'eA 

Taking the expectation and ignoring all covariances one comes to the rate equation 

dy a (t) M M 

= Va(t) }^ W aa'y a '{t)- 
a'eA 



This is exactly the rate equation (2.9). Hence, our model is equivalent to a perfect integrator with exponential 
transfer function and cumulative reset. We also want to point out that the choice of an exponential transfer 
function is well justified by physiological findings ( [51 [213] ) . 

3.4. Maximum likelihood estimation of parameters. We have already observed in the introduction that 
our models have many common features with a class of generalised linear model (T5J HH]> but see also [5TJ 122] . 
These relations are already clear from Equation (2.51. In fact, taking the logarithm of both hand-sides leads to 
the relation 

log A Q (t) := log A a (0) + N a' (t ~ c) log Waa'- 

a'eA 

This shows that the natural logarithm of the instantaneous rate is linear in the model parameters. Since the 
former is the canonical parameter of the likelihood, the relation to the GLM models is clear. We want to 
illustrate this fact with a simple computation. Assume that we want to estimate the parameters W21, W22, A, r(s) 
of a stochastic perfect integrator given the set of observations X(t), Y(t). Then the first attempt is to maximize 
the likelihood 

P[A,Y> 21 ,w; 22 ,A,r(0)] = 

[](1 -cxp(-eA)) x « Y[cxp(-e\) 1 - x ^x 

ten tew 

- exp(-er(<))) y W J[ cxp(-er(i)) 1 - y W. 
ten teu 
The input rate A does not change with time and so this is equivalent to maximize 

- exp(-er(<))) y W exp(-er(t)) 1 - y W. 

ten 
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Multiplying by pv , where N is the total number of spikes, does not change the extremal points. Moreover, one 
only has to multiply if the exponent is different from 1. All in all, after applying the usual exponential identity 
we have to maximize 

JJrfe) Yl exp(-er(i)) = JJrfo) exp(- £ er(t)) 

3 M*a 3 t^tj 

i-T 



TTr(i,)exp(- f r(s)ds). 
Jo 



3 

Applying the logarithm to both sides we finally come to the problem of maximing the expression 

(3.4) / log(r(s))dN(s) - [ r(s)ds. 

Jo Jo 

This is the natural form of the maximization condition of the generalised linear models mentioned before. This 
simple computation has two consequences: 

• the multiplicative interaction rule leads to an implementation of the generalised linear model which is 
local in time; 



• estimating the connection strengths of a neural population by means of (3.4) implicitly assumes that 



synaptic interactions have a multiplicative effect on the instantaneous firing rate. 

4. Network stability 



As we have already pointed out, numerical simulations suggest that the fixed points of Equation ( 2.9 1 correctly 
predict the asymptotic firing rate for the stochastic perfect integrator, both in the adiabatic and in the transient 
regime. Some of the data are reported in this paper, and the code is available on request. 

In the case of the SPI, the analysis was simplified by the fact that the rate equation has a single non-trivial 
stationary point. For general networks consisting of n neurons we have 2™ — 1 non-trivial stationary points; it is 
thus not immediately clear which of them arc candidates as asymptotic firing rates. It turns out that the possible 
firing rates are the ones corresponding to locally attractive fixed points of the rate equation. This finding was 
supported both from heuristic arguments and by the analysis of the activity of several different networks of sizes 
up to a few hundred neuron, both of random and engineered type. 

The natural question which arises at this point is what happens in the case of networks having more than one 
stable fixed point. In this case, the asymptotic firing rate will converge to either of the stable vectors, and the 
decision will depend in part of the initial condition and will be in part random. 

We explain this phenomenon with two examples. First, consider a network with a rate equation which possesses 
only two positive fixed points, say r\,ri- Assume further that r± is globally attractive and that r 2 is unstable. 
Then, for large times, the average activity of the network will be exactly r\, even if the network is started at the 
point r 2 in the adiabatic regime. If, instead, the network has a rate equation which possesses two positive, locally 
attractive fixed points, again r l7 r 2 , then the asymptotic firing rates of the individual neurons will be either the 
components of r\ or the ones of r 2 , and the decision will be random. We stress that this is a collective behaviour 
and that the individual asymptotic firing rates will be given by the individual components of either fixed vector, 
all components being chosen from the same fixed vector. It is not possible to observe some individual firing rates 
from the vector r% and some others from the vector r 2 . 

In this section we explore the possibility of using it to construct networks which solve certain computational 
tasks. 

4.1. General properties of the rate equation. Before we start the exploration of the possibilities of our 



model, we want to discuss some general properties of the Equation (2.9 1 
As a first step, we split our units a e A in different populations. 

Definition 4.1. We use the following notation: 
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• The set of units a for which i aa > = for all a' G A is called the input population. 

• The set of units in the input population for which £ aa = is called the pure Poisson input. 

• The set of units in the input population for which £ aa ^ is called the transient input. 

• The set of units a not in the input population form the recurrent population. 

Units belonging to the transient input can only show two different behaviours. Their activity either converges 
to 0, or explodes exponentially. For this reason we impose the following. 

Assumption 4.2. The system under consideration does not possess transient input. 

Moreover, we assume that all non-recurrent units in the system have self-inhibition. 

Assumption 4.3. If a G A is not part of the pure Poisson input, then £ aa < 0. 



Let us make an additional check for the correctness of Equation (2.9 1. Since the rate of a point process is a 
positive function, one should expect that the positive cone of K' A ' is invariant for the Equation (2.9 1. To see 
that this holds, observe that the boundary dC of the positive cone C is given by 

dC = \J{xe R lAl : x > 0, x a = 0}, 

and so invariance holds if and only if 

dy a (t) 



dt 



>0, 



whenever y a (t) = 0, but this is clear since dy ^*- ) = 0. We have just proved the following result. 



Proposition 4.4. The positive cone o/IR'" 4 ' is invariant under the flow induced by (2.9 1 



The same result holds if one substitutes the positive cone with any quadrant of the space K'" 4 ', but this is of 
course not relevant for probabilistic applications. 



As a second step, we rewrite of Equation (2.9 ) by separating the Poisson input from the rest of the population. 



To this end we denote by P d A the Poisson input of the system, define i p :— y p (0) for all p G P, and denote by 



R the recurrent population. This makes sense because y p (t) is constant for all p G P. Equation (2.9 1 can thus 
be rewritten as 



(4.1) = y r (t) I VsiWrs + 

,seR pep 



complemented with the initial condition y(0) = y$. We define Cr as the principal minor of C associated with 
R C A and Cp as the restriction of C to P C A. 

We assume during the rest of the section that the coupling matrix Cr is negative definite. In this case, it is 
in particular invertible with inverse In order for the right hand side 




} J trpip | = yr{C R y r + Cpi r ) 



of Equation (4.1 1 to vanish, we either have y r — or y r = —(C R Cpi) r . We thus obtain the following result. 



Proposition 4.5. If Cp is invertible, then Equation (4.1l has 2' A ' critical points 



Of course, not all stationary points are positive. In fact, the negative definiteness of a matrix has no implica- 
tions for the negativity of its inverse. Therefore, even in the case of purely excitatory Poisson input, it is difficult 
to draw any conclusion about the existence and number of positive critical points. 
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Example 4.6. Consider the matrix A = ( ^ ^1^' ^hen ^ ^ s negative definite, but A^ 1 = f ^ 1S 

neither positive nor negative. Assume now that the input is positive, i.e. purely excitatory. As a consequence, 
depending on the input level, each of the 2' A ' of the stationary points will, or will not, be in the positive cone. 

On the other hand, A = ( j ^l^j 1S ne S a ti ve definite and its inverse A^ 1 = ( J ^1^) 1S & ne S a ^^ ve 

matrix. In this case, for purely excitatory input all 2> ' stationary points will be in the positive cone, irrespective 
of the input level. 



Precise statements for quadratic systems like those given in Equation (2.9 1 are very difficult, see [TT] for a 
review of some open problems. However, in our case it is not difficult to see that all relevant solutions are 
bounded. To see this, define the energy function z{t) := J^a^AVa^)- An easy algebraic manipulation yields 

dz(t) 

= £r'U -y + e-y, 



dt 

for an appropriate vector e. So, ||y|| — ► oo implies dz ffl — > — oo and z(t) > because of the invariance of the 



positive quadrant for the equation (2.9 1. Summing up, if z(t) — ► oo, then dz Jp — ► — oo, and so z(t) is bounded, 



since it is positive. This proves that all y a are bounded, which is the following. 



Proposition 4.7. Assume that Lr is negative definite. Then all positive solutions of (2.9 1 are bounded. 



Although negative definitcness guarantees that solutions are bounded, the system is not dissipative. To see 



why this is the case, denote by F the right hand-side of Equation (4.1 1 and observe that 

dF r (y) 



^ £rpip ~t~ ^ ^rsVs £r 



dy r 

d pGP s£R 



■Vv 



Summing up with respect to r, we obtain 

div F(y) = ^ ( -rpip + £r.V ^rrVr- 
peP.reR reR, 

We call the three terms the (total) input, dissipation and inhibition, respectively. Of course, since Cr is negative 
(semi)-definite, one obtains the estimate 

div.F(y) < input. 

Since the dissipation and inhibition are homogeneous polynomials in y, it is not possible to replace the input by 
a better constant. Equality holds if and only if y = 0. Concluding, if the total input is positive, the system is 
neither dissipative nor conservative, although it has bounded orbits. 

4.2. General properties of two-dimensional models. We now study the simplest possible case: networks 
consisting of two neurons, each of them receiving input from a Poisson process. We assume that P = {1,2}, 
R={3,4} 

t/=(r), *33=*44 = -l. 



We further assume that the parameters £31 , £42 represent equivalent inputs and that each input unit of the input 
population is projecting to a single recurrent unit. In symbols 

l/x(*) = !te(t) = l, 4i = 4 2 = 0. 
We are analysing the ordinary differential system 

dt U/4 J U/4(-J/4 + hzVx + £42), 
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The Jacobi matrix of the system is given by 

J 

The stationary points are 



2/3 \ _ f-ty3 + 4-42/4 + 4l 442/3 
2/4 J \ 432/4 -2?/4 + 432/3 + 44 



o_ /0\ o_ /4i\ 4 / o 



y = I o J ' y = I o J ' y -V4 2 



and finally 



J/ 



l-<?34^43 

I42+I43I3I 
\ 1 — ^34^43 / 



Observe that the expressions for 2/ 1 ,2/ 2 can be easily understood intuitively. If one unit is silent, the rate of the 
other one only depends on the input fed into the active unit. The numerators of y° is also easy to understand: 
this is simply the total weight of the paths of the full connectivity matrix C leading to the corresponding neuron. 
The denominator is not as easy to understand and requires some quantitative consideration. Before we start the 
discussion of the three different exemplary cases, we make some general observations about the Jacobian matrix. 
First, denoting by a {A) the set of the eigenvalues of a matrix A, 

a(J(y )) = {4i,4 2 }. 

This means that the stability of the trivial state depends only on the sign of the equivalent input. 
In the degenerate case, i.e. when only the first neuron is active, we have 

aiJiy 1 )) = {-4i,4i4 4 + 4 2 }. 

For the second neuron the expression for the eigenvalues is analogous. 
In the symmetric case 

|4 4 | = |43| =: 4ross, 4l = 42 =: 4iput, 

the eigenvalues of the Jacobian in the critical stationary case are given by 

/ 7V c\\ r d 4 4ross J^-input -1 



4.3. Positive feedback loop. We assume that all neurons have the same self-inhibition, i.e. £n = l 2 2- Then, 
the characteristic equation of the coupling matrix is 

x 2 - 2£ u x + 4 - £ 3i g i3 = 0. 

Solutions of these equation are both negative if and only if 

4 > 4 4 4 3 , 

in other words, the network is stable if and only if the self-inhibition is strong enough. 

Silent state. In this case, 4i,42 < 0, and so the silent state is locally attractive for all possible choices of 
parameters. Observe that if the network is unstable, i.e. for large cross excitation, one can have a situation that 
for small initial values the network converges to the silent state and for large initial values activity explodes. 
Degenerate state. These states are always negative, and so they are not relevant for the discussion. In fact, if 
one of the two neurons is silent, the other neuron is not receiving any excitatory input, and so will converge to 
the silent state. This shows that no degenerate state can be stationary. 
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Critical state. Let us consider the symmetric situation where £34 = £43- First, we have to guarantee that the 
critical state is actually positive The critical question is whether the mutual excitation £^ mss ^ s larger than the 
self-inhibition la = 1. In fact, if £ cmss > 1, then y c > 0, but the network is unstable by our initial considerations. 
Summing up, we found that the positive feedback loop 

(1) is dissipative and only possesses a reachable stationary state if the self-inhibition overcomes cross exci- 
tation; 

(2) is unstable and possess a further unstable stationary point in the opposite case 




Figure 4. Illustrative examples of an oscillator with excitatory drive and of a winner-takes-all 
network. Parameters are Across — ±0.22,£i nput = 0.18, tu = —0.1. Upper panel: Simulation of 
an oscillator with excitatory drive. Observe the initial transient due to the high initial rate; the 
input rate is 20. Lower panels: simulation of two possible realizations of an winner-takes-all 
network. Parameters are Across = — 0.22, £i np ut = 0.18, la = —0.1. Input to both neurons are 
Poisson trains with rate 10. 

4.4. Oscillator with excitatory drive. In this case 

4a = -e 3 4 > 0, £31 > 0, £42 = 0, 

such that one can consider a reduced system consisting of 3 neurons. The coupling matrix is always negative 
definite, and the network is dissipative for all choices of parameters 
Silent state. The trivial state is now stable, but not attractive. 

Degenerate state. The state y 1 is positive, so it is reachable. However, it has Jacobian eigenvalues <r(J(y 1 )) = 
{—£31, £43^31}. They have opposite signs, and so the state is unstable. 

Critical state. The critical state is positive, but the sign of the eigenvalues depends on the choice of the parameters. 

It is interesting that the network is always stable; we use this example to illustrate the fact that the stability 
properties of the rate equation are equivalent to those of the stochastic dynamics. 

In the simulation plotted in Figure [4] we used the reduced system of 3 neurons with parameters 

£31 = ln(1.25), £ 2 3 = -^32 = ln(0.8), 

and finally 

^22 = £33 = -0.1 

For this choice of parameters, the critical state is attractive, so the stochastic dynamics should converge to this 
fixed point. In order to show that the predicted asymptotic firing rate is globally attractive for the stochastic 
dynamics, the initial rate for recurrent units was fixed to 1000, whereas the input rate was fixed at 20. Subse- 
quently, we estimated the average firing rate of unit 2 and 3 in the second half of the simulation and they were 
found to be 5.8 and 13.8, in good accordance with the predicted values of 7.46 and 16.65. The discrepancy is 
mainly due to the variability in the Poisson spike train used as input; in the second half of the trial presented 
in Figure [4] for instance, it actually fired only 78 spikes instead of 100. In fact, if one uses the normalized spike 
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count of the input unit in determining the firing rates of unit 2 and 3, one obtains the corrected prediction of 
5.82 and 12.99, with an error of ~ 5%. 

4.5. Negative feedback loop with external excitation. This case deserves particular attention. Intuitively, 
negative feedback loops can be used to implement winner-takes-all mechanisms, see [4 for a discussion of the 
biological relevance. The mechanism is the following: If both inhibitory neurons compete by inhibiting each 
other, the one receiving the largest part of the input could achieve to completely suppress the opponent. This 
corresponds to the situation where the degenerate state is stable. We want to analyse whether it is possible for 
the network under consideration to operate in this regime for a large range of parameters. Observe that a very 
similar, rate based model has been used in [8 , although it was derived from somewhat different considerations. 

To facilitate the analysis, let us put ourselves in the symmetric situation. We first observe that the network 
is not dissipative, since the coupling matrix is not definite, and the input is positive. However, all solutions are 
bounded. To see this, observe that every neuron receives bounded excitatory input, and so the output rate is 
bounded by — j^Ui for the first neuron and by — f° r the second one. 

Silent state. The state is repelling since both eigenvalues of the Jacobian are strictly positive. The first conclusion 
is that the activity of this network will never fade out. 

Degenerate state. Let us start by neuron 1. One of the eigenvalues is always negative. The second is negative 
if and only if £31 > jf^y- In other words: the degenerate state of a neuron is attractive if and only if the own 
input overcomes the input of the other neuron divided by the cross inhibition. The same happens for the second 
neuron. We have two distinct regimes: 

(1) The cross inhibition is larger than the self-inhibition, i.e. ^34 1 > 1. In this case at least one of the 
degenerate states can be attractive, depending on the level of the input. In some case, both degenerate 
states can be attractive. 

(2) The cross inhibition is smaller than the self-inhibition, i.e. | ^34 1 < 1. In this case at most one of the 
degenerate states can be attractive, depending on the level of the input. This means that if the inputs 
are close, the network will not converge to a degenerate state. 

Critical state. An easy computation shows that the critical state is attractive if | ^34 1 < 1, and unstable otherwise. 
Summing up, for low levels of cross inhibition, one could have that the critical state and possibly one degenerate 
state are attractive, depending on the input level. 

Conversely, for high levels of cross inhibition, the degenerate state corresponding to the neuron receiving the 
most equivalent input is always attractive and possibly also the second degenerate state. 

This means that the actual stochastic trajectory will end up in one or the other state, depending on the 
realization. In fact, for high level of cross inhibition and for inputs which are close, the "winner" will be chosen 
randomly, and the probability depends on the actual ratio of the inputs. 

To illustrate this phenomenon we show two different simulations for the same input level with different 
outcomes. In Figure|4]one can appreciate the stochastic properties of the network. Although the initial conditions 
of the network are exactly the same, the system evolves into two different states, each of them corresponding to 
one of the attractive, degenerate states of the network dynamics. 

We also want to point out that even in the framework of minimal networks like the ones we have just 
investigated, the networks showing the most interesting dynamics are those which possess inhibitory neurons. 
Again, we stress that this possibility is not given for networks of Hawkes's processes, such that our model really 
represents an important step toward modeling and understanding the dynamics of neural networks. 

5. Discussion and outlook 

5.1. Summary. We introduced a class of stochastic point processes with multiplicative interactions, where 
dynamic changes in the rate of each component process are induced by the events in all other processes connected 
to it. 
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We chose multiplicative interactions both for biological and mathematical reasons. From the biological point 
of view, the model obtained corresponds to a (non-leaky) integrate-and-fire neuron with linear synapses and 
exponential transfer function. From the mathematical point of view, multiplicative interactions allow us to 
treat inhibition, without explicitly invoking membrane potential dynamics, but leading to a formalism that is 
analytically tractable. 

We outlined the general theory of such systems and proved that important aspects of their temporal evolution 
are described by a differential relation involving expectations, covariances and infinitesimal terms of first order. 
We could make some first steps in elucidating the relation between the deterministic dynamical system of expected 
rates and the stochastic dynamics of the interacting point processes. In fact, extensive numerical simulations of 
networks of different sizes and architectures, as well as some heuristic analytical arguments, clearly indicate that 
the stability of the stochastic point process dynamics is equivalent to the local attractiveness of the corresponding 
fixed point of the rate equation. 

We then moved on to the analysis of a stochastic perfect integrator model. We illustrated to which degree the 
fixed point of the rate equation predicts the firing of a stochastic point process in equilibrium. A master equation 
for the time evolution of the rate distribution was derived, and we supported our findings by simulations and 
numerical results. The master equation shed some new light on the equilibrium distribution of the rates. In 
the future, it could also be used to extract information about the transients, but so far we did not attempt to 
actually find solutions to the equation. 

We compared our multiplicative model with related approaches in the neuroscientific literature, formally 
proving that 

(1) our model corresponds to an integrate-and-fire neuron with linear synapses and exponential transfer 
function, 

(2) it is a generative model for the framework described in [T5] , 

Finally, we analysed the differential system of the rates in some simple, biologically relevant cases. It turned 
out that it is possible to easily implement a robust winner-takes-all decision mechanism for this type of networks, 
similar to a rate model that was introduced previously based on heuristic arguments [8]. We first studied the 
stability properties of the equation analytically, and then presented simulations that confirmed the empirically 
observed equivalence of stochastic and deterministic stability of fixed points. 



5.2. Adiabatic and transient regimes. We have already pointed out that the rate equation appears to 
correctly predict the behaviour of the stochastic system only in the adiabatic regime. However, this concept 
is not completely specified and must be investigated further. We see three different possible approaches to the 
problem: 

(1) deriving a rate equation for the covariances of the rates based on the master equation, and showing that 
they all asymptotically vanish; 

(2) developing a quasi-Floquet theory for stochastic systems and deriving conditions under which a trajectory 
of the stochastic system converges to a trajectory of the deterministic system; 

(3) employing abstract Martingale theory to develop a genuine probabilistic approach to interacting point 
process dynamics. 

All three approaches are mathematically challenging, and it is not clear whether they can be successful given 
the current state of mathematical techniques. 

Even more challenging is the issue of transient behavior. Stochastic transients are highly relevant for neural 



signal processing, but mathematically difficult to analyse. In Section 3.1 we have seen how the immediate 
response to a step input, given a constant initial rate, exceeds the equilibrium response. Such a mechanism could 
contribute to phenomena like population spikes in the auditory pathway. Transients are of course not specific to 
our model, but common for many non-Poissonian point processes, e.g. for renewal processes with positive ageing, 
if the time-dependence of the hazard rate is arranged properly [H] , 



INTERACTING POINT PROCESSES 



17 



In principle, it is possible to understand the transient behavior of the stochastic perfect integrator via its 
associated master equation. However, this equation is difficult to solve analytically, and this makes it difficult to 
extract information from it. 



5.3. Specific neuronal circuits. The computations explained in Section |4~2| have shown that a negative feed- 
back loop can be used to set up winner-takes-all networks, with excellent performance in the low-rate regime. 
It is an interesting question, whether the performance of such circuits can be related to specific parameters of 
the corresponding deterministic system. Possible candidates are the Lyapunov exponents, or some measure for 
the size of the basin of attraction. A related question is how to induce alternating behavior, as observed in 
common models of binocular rivalry, see [7] . Preliminary investigations have confirmed that it should be possible 
to construct a model of competing neural populations that relies on the same architecture as the one described 



in Section 4.2 which reproduces many characteristic phenomena, and which is analytically tractable. 

Another important issue in the context of specific circuits is the role of global inhibition for the stabilisation 
an excitatory network. We have already seen in Section |4.2| that an oscillator with excitatory drive is always 
stable, independent of the parameters. Further simulations (not shown) suggest that global inhibition has a very 
good stabilizing effect on excitatory networks that are otherwise unstable. A study of this problem reduces to 
the spectral estimation for the special type of matrices corresponding to the circuit under consideration. 

5.4. Random networks. A study of large random networks should be performed; as a matter of fact, a crucial 
test for the model is whether it is able to reproduce statistics of parallel spike trains as observed in cortical 
recordings. The classical approach pQ is to derive a self-consistent equation for the parameters under investigation 
and solve it to characterize the states in which the network can operate. 

Important progress in this direction has recently been achieved |21j . In fact, the type of mean- field approx- 
imation worked out in that paper relies on some type of randomness in the underlying network; the authors 
derived an ODE system for the time evolution of mean rates and covariances, and they showed that it correctly 
predicts the network behavior. 

The advantage of our approach is that it is possible to explicitly include the topology of the underlying 
network into the description of its activity. Taking into account the issues that we have discussed in Section |5.3[ 
we want to explore the possibility of embedding specific neuronal circuits into some appropriate class of random 
networks. The final goal would be to understand the computations which can be performed by biologically 
structured random networks. 

5.5. Extensions of the model. Our model can be extended into different directions. First, reasoning as in 
Section [3. 3[ one could derive rate equations also in the case of leaky integrate-and-fire neurons, or one could add 
a refractory period to the single-neuron dynamics. Preliminary studies in this direction have been performed, 
which show that rich behavior arises, including periodic trajectories of the population activity. 

Obtaining information about the time evolution of higher moments is also of great importance to correctly 
address the issue of transient behaviour. In Section |3.1| we have shown how to derive a system of differential 
equations for the moments and, if a master equation for networks can be derived, the same method could be 
employed to derive a system for higher moments of networks. 

Acknowledgments. We wish to thank an anonymous reviewer of the paper for discovering a flaw in our derivation 
of the master equation, and for suggesting the appropriate correction. This work has been supported by the 
German Federal Ministry of Education and Research (BMBF grant 01GQ0420 to the BCCN Freiburg). 

Appendix A. Elementary facts about infinitesimal random variables 
We start computing the conditional expectation. By definition 

E[X | r] = 1(1 - exp(-re)) + 0(exp(-re)). 
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According to the definition of a derivative, 



E[X | r] = 1 - exp(-re) 

= — (exp(— re) — exp(O)) 

dx 

= er + o(e) 



■ exp(— rx) 



e + o(e) 



z=0 



We conclude that 



EX = ^ E[X | A]P[r = A] 

A>0 

= ^ e(A + o(e))P[r = A] = eEr + o(e 2 ). 



A>0 



To see that also formula (2.2 1 holds, observe that 

exp(re) = 1 + re + o(e), rel 

So, it is apparent that 

E(l - exp(-re)) = eEr + o(e). 



Formula (2.3 1 follows from 



Var(X) = EX 2 - E 2 X = EX — E 2 X = E.Vl l - EX). 



Appendix B. Expectation relations 



Derivation of formula \ 2. 7\ We fix an arbitrary time t and compute by the formula (2.6 1 

A\ a (t) _ AA a (0)exp(£ a , eA A a ,(t-e)loguw) 



At 



A a (0) 



At 

Al\ a , eA exp(N a ,(t - e)logw aa ,) 
At 



Using now (2.4 1, the latter equals 



X i(t) 

A o (0) J| exp(N a ,(t - e)logw aa ,) ^ — — logw aa 



a'eA 



a'eA 



= K{t) Xa — W aa > ■ 

a'eA 



Because of relation (2.1 1 



E Y, ^^loguw = ^a>(t)logw aa , +0(e). 



a'eA 



a'eA 
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Computing 



E 



AXgjt) 

At 



A„(t) 



E 



A„(t) ^^log«w | A B (t) 



a'eA 



= A,(t)E 



E ^^logiiW | X a (t) 
i'eA 

X a (t) I EXa'(t)logw aa ,+0(e)\ 

Va'eA / 



completes the proof. 



□ 



Derivation of formula \2.8\ For each path of the stochastic process the fundamental theorem of calculus implies 

AA a (s) 



So, by linearity of the integral 



A„(t) = A(0) 



EX a (t) = A(0) 



[0,t]m 



At 



-ds. 



[o,t] H 



As 



Interpreting EA a (f ) as a function of time, the above relation means, again by the fundamental theorem of calculus 

AEA Q (i) _ E _AA a (t) 



At 



At 



We compute as in the first part of the proof of Equation 2.7 to obtain 



E ^ = E 

At 



Aa(i) E 



logw Q 



a'eA 



By the linearity of the expectation, the latter satisfies 

*«'(*) 



E 



Aa(t) E 



log 1U Q 



a'eA 



: ^ E L(f)^Wl0 gU , 00 , 

a'GA L 

: E log«w(E[A a (t)A„,(f)] 

a'eA 

•0(e)EA (t)). 



We have to justify the last equality. First, 



E[A a ^^] = E p [ A -' = M]E[Aa^|A Q ' - l4- 
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By the conditional independence of X a i and A a , and by (2.1 1, the latter can be written as 



M>0 



H>0 



£>[A a , =/z]E[A a |A a , = M ] (M + 0(e)) 



Ai>0 



: 53 ( M + 0(e)) i/P[A a / - /x]P[A = i/|A„/ = /i] 
: E[A a A a /] + 0(e)E[A Q ] 



□ 



Appendix C. Derivation of the master equation 

Let us denote by B(r,6) a ball centered on r and of radius S. Because of the conditional independence of 
X(t),X(t + e),Y(t),Y(t + e) we obtain that 

P[r(t + e) G B(r,S)] = 

r 8 
W21W22 ' W21W22 
r 8 



r(t) G B 
r(t) G B 



r(t) G B 



W21 w 2 i 
r 8 \ 



P[X(t) = l)V[Y(t) = 1] 

¥[x(t) = i]P[y(t) = 0] 



p[x(t) = o]p[y(t) = 1] 



\W 2 2 W22/. 

r(t) G B(r, S)]P[X(t) = 0]P[K(t) = 0] 
1 



W2lW 22 



r(t) G B 



W21W22 



x (1 — exp(— Ae))(l — exp(- 



W21W22 ' 



-)) 



W22 



r(i) G B I — ,5 
w 2 i 



(1 — exp(— Ae)) exp( e) 

W21 



r{t) G B 



W22 



,<5 



exp(— Ae)(l — exp( e)) 

W22 



\r{t) = r] exp(— Ae) exp(— re) 



The last term can be written as 



So, defining 



P[r(t) = r] + P[r(i) = r](exp(-Ae) cxp(-re) - 1). 
f{r,t) :=P[r(t)6 5(r,i)] 
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for a linear infinitesimal 5, and rearranging appropriately, we come to the relation 



df(r, t) 1 r 1 — exp(— Ae) 



dt 



"/(- 



-(1 - exp( 



W21W22 W21W22 

1 , r 1 - cxp(-Ae) r 

H / ,*) exp e) 

W21 W21 e W21 



W21W22 



1 r 1 - exp(- 

+ — / — ,t)exp(-Ae — 

W22 W22 e 



:<0 



/(r,t) 



exp(— Ae) exp(— re) — 1 



We apply now the usual exponential identity, and ignore all infinitesimal terms to come to the differential equation 

^ = —/(—,*) + ^2-H — ,*) - (A + r)/(r,t), 

OT W 2 1 W21 W22 W22 



which is exactly Equation (3.1l. 



Appendix D. Derivation of the moment equation 

The moment equation has been derived following the suggestions of an anonymous reviewer of the manuscript. 
To see how it works, remind that the n-th moment is defined as 



So, deriving with respect to time 



dt ~ dt 



/*»(*) := / r n f(r,t)dr. 
Jo 

r n f(r,t)dr 
,df(r,t) 



dt 



-dr 



r n — f( — ,t)dr+ r r n ^f( — ,t)dr 

W21 W 2 1 Jo W22 W22 



OO 



-X r n f(r,t)- r n +'f(r,t)dr 
Jo Jo 

poo \ 

= w 2 i / s n w^ 1 f(s,t)ds 

Jo w 2 i 

+ W22 / s n w% 2 —f(s,t)dr 

Jo w 22 
- \(l n (t) - (Jtn+l(t) 

= \w 2 \^n{t) + W22^n+l(t) - A/i n (t) ~ fX n+ l(t), 

where we have repeatedly performed integration by substitution. A more compact writing is 

d ^ = K 2 - l] Mn+ i - A[l - w^ n (t). 



This is Equation (3.3). 
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